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Abstract 

An implicit finite-difference algorithm is developed for the num- 
erical solution of the incompressible three-dimensional Navier-Stokes 
equations in the non-conservative primitive-variable formulation. The 
algorithm is second-order-space-time accurate, iterative and enploys 
a point SOR method for solution. 

The flow field about an airfoil spanning a wind-tunnel is computed. 
The coordinate system is generated by an extension of the two-dimension- 
al body-fitted coordinate generation techniques of Thompson, as well as 
that of Sorenson, into three dimensions. Two-dimensional grids are 
stacked along a spanwise coordinate defined by a simple analytical 
function. 

A Poisson pressure equation for advancing the pressure in time is 
arrived at by performing a divergence operation on the momentum equa- 
tions. The pressure at each time-step is calculated on the assumption 
that continuity be unconditionally satisfied. An eddy viscosity coeff- 
icient, computed according to the algebraic turbulence formulation of 
Baldwin and Lomax, simulates the effects of turbulen'~e. 
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Computational Raaulta ara praa«itad for tha KACA 0012 and M A CA 
66^018 airfoils in sloulatad wlnd-tunnal tasts. Computations vara 
parformed for a Reynolds numbar of 1000 for tha HACA 0012 airfoil at 
angles of incidence of 0 and 12 degrees. They Reynolds number was 
increased to 10,000 for computations with the HACA 66^018 airfoil as 
the profile under investigation. Pressure plots and velocity profiles 
are presented for the various cases. Comparison with wind-tunnel data 
for the HACA 663 OI 8 airfoU indicates qualitative agreement with typ- 
ical pressure distributions. Quantitative replication of experimental 
data was beyond the scope of this dissertation, what with scant data 
points in the discrete mesh systoa and the consequent inability to 
perform computations at high Re 3 ntolds numbers. The size of the grids 
and the total computation time were restricted by considerations of 
affordability. The feasibility of an implicit solution of the complete 
three-dimensional Havier-Stokes equations is established and sufficient 
encouragement for future computations at realistically high Reynolds 
numbers is derived. 
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Chapter 1 


Introduction and Literature Review. 

The iiiq>ortance of and the necessity for conq>uter solutions of 
the complete three-dimensional Navier-Stokes equations as a reliable 
alternative to experimental fluid-dynamic studies cannot be over empha- 
sized. It was, therefore, not the lack of intellectual motivation and 
initiative but rather that of sufficient advances in computer technology 
that limited the researchers' involvement in this field since the 
beginning of the computational fluid dynamic era in the early 1960's. 
Graves [l] has given a comprehensive overview of the evolution of the 
CFD discipline. Until the advent of the supercomputers the CPU time 
as well as storage requirements associated with three-dimensional 
Navier - Stokes solutions were truly prohibitive. These difficulties 
have been partially overcome by recent computers with radically changed 
architectures permitting enormously Increased speed and storage capacity , 
particularly the CDC 200 series and the CRAY 1 machines. 

The two chief tasks comprising fluid dynamic computations, namely 
effective grid generation and developing accurate algorithms for solving 
the governing equations play complementary roles to each other. Con- 
sequently, advances in the techniques in these two fields have gone 
hand in hand. The concept of boundary-fitted coordinate systems developed 
by Thompson, et. al. [2], has established itself as particularly suit- 
able for effective grid generation for arbitrary shaded bodies. 

Sorenson , et. al. 13] , offer a different version of the body-fitted 
grid generation technique. Both these techniques are extensible into 
three dimensions for bodies with simple geometries without necessitating 
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any rigorous mathematics. Boundary* fitted coordinate systems permit 
computation in a rectangular grid with straight boundaries in the trans* 
formed plane, thus making the in^)lementatlon of boundary conditions 
much simpler than in the physical plane with arbitrarily shaped 
boundaries • 

A primitive variable formulation of the complete Navier*Stokes 
equations is advantageous for three-dlmr.. clonal flow computations. 
Various methods, each with its own area of applicability, have been 
formulated and successfully executed by researchers in flow situations 
of different characteristics. These methods can be broadly categorized 
into explicit and implicit algorithms; their relative advantages being 
less computer storage requirements and greater speed respectively. 

Explicit techniques are restricted by time-step size limitations 
dictated by stability considerations. Among implicit techniques the 
Approximate Factorization technique, formulated by Beam and Warning [4J, 
has been successfully used in computation of compressible flow, but 
is not readily applicable to incompressible flow. However, modified 
versions of AF techniques have been used by Steger and Kutler [5J and 
Bernard [6j for attempting incompressible flow computations. Given 
the above considerations and the task of solving the conq)lete three- 
dimensional Navier-Stokes equations, a fully implicit scheme employing 
successive over- relaxation (SOR) emerged as the logical alternative. 
Available information on numerical solution of three-dimensional 
incompressible turbulent flow is scarce. East and Pierce [7], 
provide one of the very few publications in this field* 

The three-dimensional coordinate systems for the problem at hand 
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are generated by means of extensions of the 2-D boundary-fitted grid 
generation techniques of Thoiiq>son [8] and Sorenson [3]. The spanwiae 
extent of the coordinate syston is achieved by stacking a nuad>er of 
2-D C-type grids along the span of the wing. Coordinate lines are 
made to concentrate in the boundary layers on body surfaces, i.e., the 
wind tunnel walls as well as the airfoil. Computatims for the veloci- 
ties and the pressure are performed in the transformed plane. 

the governing equatiMis are the incompressible Navier-Stokes 
equations written in their non-conservative forms in terms of the 
primitive variables. A Poisson equation for pressure is obtained by 
taking the divergence of the momentum equations. The eddy viscosity 
in the Reynolds -averaged momentum equations is derived from a two-layer 
algebraic turbulence model as proposed by Baldwin and Lomax (9]. 

Boundary conditions imposed on the walls of the wind tunnel and the 
airfoil are obtained employing no-slip conditions for velocities and 
the normal derivative of pressure from the momentius equations evaluated 
on the boundaries. Constant free-stream conditions are Imposed on the 
upstream boundary. The flow is initiated by using a cosine body force 
to gradually accelerate the airfoil together with its attached coordi- 
nate system from rest to the free-stream velocity over a period of time. 

A fully implicit algorithm is obtained by representing the non- 
dimensional Navier-Stokes equations by second-order backward-time and 
central-space finite difference approximations in the transformed 
coordinate plane. The resulting coupled system of nonlinear algebraic 
equations is solved by means of a point successive over-relaxation 
(SOR) iterative method with locally optimum acceleration parameters 
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The development of the algorithm was along the lines followed by 
Thomspon [13] in two-dimensional computations. 

The succeeding chapters of this dissertation deal in further 
detail with the coordinate generation techniques, the formulation of 
the governing equations and the implementation of the numerical tech- 
niques in that order. It was the desire of the author to present the 
essence of the elements Involved without burdening the reader with 
cumbersome repetitions. Solving the ccmplete three-dimensional Navier- 
Stokes equations should be tantamount to generating and analyzing all 
the flow phenomena arising out of a given fluid-body configuration. 

It is the hope of progress towards achieving this cherished goal of 
researchers in computational fluid dynamics that motivated this research. 


Chapter II 


The Coordinate System. 

2.1 The boundary-fitted concent. 

A proper coordinate system can be a quite unforgiving prerequisite 
for the ease and accuracy of solution of a physical problem. The pro- 
priety consists mainly of freedon as to the shape of the body being 
investigated. Realistic aerodynamic problems rarely involve body sur- 
faces amenable to analytical representation and often Include sharp 
comers implying discontinuous slopes. A Cartesian coordinate system 
under such circumstances would be limited in efficacy by the need for 
Interpolation near the boundary surfaces. 

The boundary-fitted coordinate system of Thompson, et. al. [8], 
provides a set of curvilinear coordinates with one line of each family 
coincident with a boundary in the physical domain. This eliminates 
the need for interpolation and results in a complete freedom regarding 
the shape of the bodies involved. It has the added capability for 
concentrating coordinate lines in regions where large gradients in the 
physical, quantities are expected to arise. 

2.2 Theoretical Formulation 

Computational grid generation in its essence consists of a mapping 
between real and computational spaces. The physical space defined by 
Cartesian coordinates x and y is mapped onto the computational space 
through the mapping functions 

5 - C (x, y) (2.1) 

and n ■ n (x, y) 

5 



One-to-one correspondence is ensured by subjecting ^ and n to the suff- 
icient, though not necessary, conditions of satisfying the extremum 


principle and of monotonic variation in 5 5. ^min — — 

n . The extremum principle precludes occurrence of extrema in ^ and 
max 

n within the field and is complied with by making the inner and the outer 
boundaries coincide with the n ” and the n " lines respective- 

ly, while £ . and C occur on the outflow boundary, 
min max 

The topological correspondence in a C-type grid about a 2-D air- 
foil may be better understood with the help of Figure 1. The boundary 


^ " ^max is mapped onto the inner boundary ^'5 “ containing the 

branch cuts and the airfoil. 5 ■ 5 and C ■ £ correspond to the 
downstream sections and respectively, 5 increasing clockwise 
around the airfoil. The n ■ constant family of lines form open curves 
resembling the letter C. The n - boundary is mapped onto the 

outer boundary T2 - - F^ comprised of the wind tunnel walls and the 

inflow section F^ 

With the afore-mentioned desired characteristics of ^ and n in 
mind, it is natural to envision them as solutions of elliptic partial 
differential equations with Dirlchlet boundary conditions. Elliptic 
equations permit any desired distribution of C and n on the boundaries. 
Moreover, the Inherent smoothness of solutions of elliptic systems is 
well recognized. The choice of elliptic systems is further reinforced 
by the ability of the Inhomogenous terms in Poisson's equations to 
control coordinate line spacing with respect to curve or a point with- 
in the field. The generating system of Poisson's equations has the 


following form: 
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5 + C 

XX yy 


P(C,n) 


\v “ T Q(5,n) (2.2) 

XX yy 

c 

A judicious choice of the distributions of P and Q makes it possible 
to concentrate lines in any desired region. Section 2.3 deals with 
actual forms of P and Q in further detail. An interchange of dependent 
and independent variables enables one to perform all computations in the 
transformed field. The generating systaa in the transformed field be- 
comes : 

“c*55 ■ ''c*nn " > 

Vii - ^*chn * ''c>’nn ' 

the boundary conditions being specified for x and y by the known shape 
of the body surfaces. The coefficients in the system 2.3 have forms 
given below: 

a ■ x^ + y^ (2.4) 

c n n 


■ x_x + y.y (2.5) 

^c 5 n •'5'n 

■ »{ + y\ (2-6) 


■ Vn ■ Vn 

The increased complexity of the transformed equations (2.3) compared 

with equations (2.2) is far overshadowed by the ease of Implementing 

boundary conditions on straight boundaries in the transformed plane. 

A great deal of simplification in computation results if integer 

values are assigned to C ^d n. To this end, increments AC and An are 
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chosen to be unity by construction. This gives rise to uniform specing 
in the transformed plane where the values of AC and An » rather than 
those of C end n are relevant to the finite difference expressions. 

The generating system of equations (2.3) is represented by second>order 
finite difference approximations in the transformed plane. The quanti> 
ties AC and An disappear by cancellation in all difference equations. 
The resulting systan of non-linear difference equations is solved by 
means of a point SOR iterative scheme. 

2.3 Coordinate Line Concentration 

It is Imperative to concentrate coordinate lines in regions with 
large gradients in flow variables, e.g., near body surfaces. Thompson 
et . ai. [10], and Sorenson achieve this effect in two different ways. 


2.3.1 Thompson et. al. approach 

Thompson's approach consists of finding a correspondence between 
n values and the radii of concentric circles distributed between two 
circles with radii r^ and r 2 > one circumscribing the airfoil and the 
other tangential to the outer boundary respectively. Applying the 
coordinate generating equations (2.3) to the radii r^ and r 2 » while 
noting that n ■ 1 on the airfoil and n ■ J on the outer boundary, 
results in the following expression for Q: 


where 


and 


II 1 

Q(C.n)- - (-, - ^) 
r 


r(n) • r 2 


+ G(n) - 1 
G(n) + 1 



G(n) 


b + r- 


^b - r. 



( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 
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The effect of Q, thus defined, is to place a line corresponding to 
n ■ K at a distance proportional to r^^ - r^^ from the body surface. 

In order to ensure proper resolution of the boundary layer the first 
line away from the boundary is placed at an approximate distance of 
one percent of the Blasius flat-plate boundary layer thickness from 
the body, i.e., 

'(n - 2) - 'l ■ 

2.3.2 Sorenson *8 Approach 

Sorenson [3j determines the control functions P and Q iteratively 
to produce specified spacing S of the first ri**line from the boundary 
and a specified angle of Intersection 6 of C-llnes with the boundary 
(fig. 2). 

P and Q are defined as: 


P(5, n) ■ p(e)e"®’^ + r( 5 )e’‘^^’’max 
Q(5, n) - + s(5)e*‘*^'’max ' 


(2.11a) 

(2.11b) 


where a, b, c and d are positive constants. The coefficients of r 
and 8 become vanishingly small at the body (n ■ 0) , so that 

P(e. 0) - p(0 


Q(C, 0) • q(5) 

Similarly at the outer boundary (n • > 

Vx> • 


( 2 . 12 ) 




(2.13) 
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Combining equations 2.12 with the generating system of equations 2.3 
we get 


P(0 




1 


q(0 


Jc 


1 (2.U) 


where, 


-(oX,. - 2BX. + oX ) 

f Si SD_Q!1 1 


«2 


-(oY,. - 2BY,. + Y ) 

,[ SS 5lL_niL 1 

j 2 ^ 

c 


n ■! 


(2.15) 


Similar expressions for r(^) and s(C) are obtained from equations 2.13 
and 2.3. 

The steps Involved in the iterative method for generating the 
grid can be summarized as follows : 

1. The four geometric constraints of specified spacing ai.d angle of 

intersection at the inner and outer boimdarles translate into equations 

which can be solved for all the derivatives occurring in equations 2.14 

except X ^ and Y on the boundaries, 

nn nn 

2. Assumed initial values of p, q, r and s are used to determine P 
and Q in Che field. 


3. The elliptic generating system, (2.3), is solved for X and Y in 
the field. 


4. X and Y are computed on the boundaries. 

nn nn 

P> q> 8 snd subsequently P and Q are determined from equations 
2.14 and 2.15. 

Steps 3 to 5 are Iterated to convergence. 
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2.4 The 3«D extualon. 


So far the generation of 2-D grlda haa been diacusaed. The body 
under investigation is a 2-D wing spanning a wind tunnel and therefore! 
an extension of the 2-D grid into three dimensions is necessary. This 
has been achieved in a straight forward way by placing 2-D grids at 
different stations along the span as Illustrated in figure 3. Since 
no change in airfoil geometry occurs along the span, it is sufficient 
to perform flow calculations on only one side of the mid-span section. 

An extra grid has been placed beyond the mid-span section for storing 
variables in order to facilitate finite differencing. 

The 2-D grids are identical to each other so that X and Y are 
independent of Z, the Cartesian coordinate along the span. Consequently, 
Z is a simple analytical function of C, 

Z - F(C) 

C, being the third mapping function Introduced by the 3-D extension. 

F can be chosen to be an exponential function in order to effect 
a desired distribution of the grids along the span. 
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Chapter III 


The Governing Equetlone 
3.1 The Nevler-Stokee equetlone 

The Nevler-Stokee equations, derived froa Newton *e law of conserv- 
ation of moaentua end Stokes' lew regarding relationships between stress 
and strain In a eater ial fluid voluae, remain to date the most versatile 
descriptor of fluid aotlon. The governing equations in this study are 
the time-dependent, Incoapresslble, Reynolds averaged Navler-Stokes 
equations in three dimensions formulated in terms of the primitive 
variables. 

In their unabridged fora the Navler-Stokes equations are as they 


appear below: 


3(pu.) 3(pu. u.) 


32 . 


3t 


ifi. + L- 

3t 3x 


iji . . |£- + :;ii + pg 

3xj 3Xj ^*1 


(3.1) 


(3.2) 


The Index 1 denotes the Cartesian diiections Xj^, X 2 and x^; and J Is 
subject to the Einstein suaaatlon convention. The quantities represent- 
ed in these equations are 



density 

velocity component 
Pressure 

shear stress tensor 


3u 3u 

^ ax, 


)-yU 


3x. 



body force per unit maes 
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(3.1) rcprescnta the nonmtuB aquatloae governing fluid flow while 

(3.2) ie the equation of continuity. Theae equetione ere direct diff- 
erential deecendente of their integral foms in teraa of conaerved 
quantities. In this "conservative" fora the equations are particularly 
effective In flow regions containing discontinuities. Since the flow 
under investigation in this research contains no such discontinuities, 

a further sinpllficatlon can be effected by incorporating the continuity 
equation (3.2) into the noaentum equation (3.1) in the following Banner. 

Expansion of all the derivatives in equation (3.1) and rearrange- 
ment of the quantities while noting tliet the density p is e constant in 
incompressible flow reduce the equations to 


at 


+ u 


j ax, 


a^ 

“i axj 


1 ^ 
p ax 


3u. 




* au, 3u, 

4- (—^ 4- — ^) 1 -f a 

^ axj ^axj ax^''^ »i 


(3.3) 


and, 


3u 

ax 


(3.4) 


J 


The third terms on both sides of the sign of equality in equation (3.3) 
vanish by virtue of equation (3.4), implying that continuity is assumed 
to be unconditionally preserved in the flow. The products of this 
manipulation are Che following "non-conservative" governing equations: 


9t j 3x. 


3u 

lx 


1 4 . 4 . . 

p ax^ “i ^ ^ 


(3.5) 


(3.6) 
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Additional cogency Is lent to the arguments for choosing the non- 
conservative form of the governing equations by the fact that nullity 
of the divergence of velocities, l.e., continuity. Is a requirement 
essential to the scheme for advancing the pressure In time as will be 
demonstrated In a latter section In this chapter. 

Finally, the governing equations are rid of dimensions by substi- 
tuting the flow variables with the following non-dimensional quantities: 



(3.7) 


(3.8) 

tu (L 
00 

(3.9) 

v/w. 

(3.10) 

(P - p )/(pU^) 

(3.11) 

CM 8 
00 

(3.12) 


where the characteristic quantities are the free-stream velocity U^, 
the airfoil chord length t, and the free-stream viscosity U^. 

The resulting non -dimens ionallzed Navler -Stokes equations governing 
inconiDressible flow are: 


~ .. 3u . - .. „ 3u. 3u 

3u ^ i 9P . 1 r „2 . , 1 _i. j\ i . 

+ u __ + -r-luV u + -77- + -:r^)j + g. 

3t ^ 9x. , e 9x. 3x. 3x. 

J i j j 1 


(3.13) 


and 


3u. 

---i =. 0 
3x^ 


(3.1A) 
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R denotes the Reynolds 
0 

number . 


The remainder of this dissertation deals exclusively with non-dl3ienslonal 
quantities and therefore the dlereses (*') will not appear henceforth. 

In usual Cartesian notation the governing equations (3.13) and 
(3.14) translate to, 

Ut + uu^ + vUy + wu^ - - + 2y^u^ + Uy(Uy + v^) 

+ + *1 (3.15) 


V..+UV + w +WV “-P + + 2pv + u(v + w) 

t X y z yR^‘*^ yy zzy 

+ Ujj(Uy + + g2 (3.16) 

w*. + uw„ + vw + ww_ ■ - P_ + w + 2p w + u (w + u_) 

t X y z zR'' zz x x z 

e 

+ My(Wy + Vjj)] + gj (3.17) 

and, 

Ujj + Vy + ■ 0 (3.18) 

3.2 Turbulence modelling . 

The effects of turbulence are simulated by appending an eddy vis- 
cosity coefficient to the molecular coefficient of viscosity p. The 
scheme for determining Is patterned after that of Baldwin and 
Lomax [9], who devised a two-layer algebraic turbulence model for sep- 
arated flo'-'s The non-dlmenslonal viscosity coefficient In equations 
(3.15) - (3.17) is replaced by 

p - 1 + e (3.18) 
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where e Is the ratio of to the molecular coefficient of viscosity. 
A few salient features of the turbulence model are: 

1) The need for finding the edge of the boundary layer is obviated by 
using the distribution of vortlclty to determine length scales. 

2) An Intermit tency factor Is employed In order to modify the eddy 
viscosity coefficient In transitional regions. 

3) Transition Is assumed to occur at the points of minimum pressure 
on the upper and the lower surfaces of the airfoil. 


3.3 The pressure equation. 

A direct temporal discretization of pressure Is not possible be- 
cause of the absence of a time-derivative of pressure in the governing 
equations (3.15) - (3.17). Advancement of pressure In time, unlike that 
of the velocities, is therefore not a straightforward task. The 
difficulty can, however, be circumvented by a stratagem of determining 
the Instantaneous pressure subject to the constraint of satisfying 
continuity throughout the flow field. In other words, the pressure 
adjusts Itself In order to force the divergence of the velocity vector 
to zero. Implying a strict conservation of mass. 

To achieve this purpose a divergence operation is performed on the 
momentum equations (3.15) - (3.17) to arrive at the following Poisson 
equation for pressure. 

2 2 2 2 2 2 2 
D^ + VP»-(u + v +w +2vu + 2wu +2wv) + — [y 7 u + y 7 v 
t xyz xy xz yzR^x y 

2 

+ y7w+y y + y v + y w + y (u + v)+y (w + u) 
z XX X yy y zz z xy y x zx x z 


+ y, (v + w )] 
zy z y 


(3.19) 
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where D le the divergence, **]( '*' *z* velocity vector. 

Several variations on this theme have been presented by Hodge [11] 
and Shanks [12]. In the non-conservative equation (3.18) is expected 
to retain an appreciable value even though D is assumed to be analyti- 
cally zero. From a computational point of view serves as a correct- 
ive term that adjusts the pressure in an effort to enforce the contin- 
uity constraint. 

3.4 Coordinate transformation. 

The governixig equations must be transposed to the computational 
grid before performing the finite difference approximations on the 
derivatives. The spatial derivatives must be transformed by the in- 
clusion of terms, relating the discrete mesh to the physical grid, in 
order to remove the physical coordinate systma from the problem. 

The derivatives, u , u , etc. in equations (3.15) - (3.17) are 
therefore substituted by their values computed according to the express- 
ions for the derivatives in the transformed plane listed in appendix A. 


Chapter IV 


The Boundary Conditions. 

The solution of the Navier-Stokes equations la a quintessential 
boundary value probles. Extroae care must be exercised in specif ing 
the boundary conditions, since the validity of the solution Is depen- 
dent on the accuracy of the boundary values. Presented below are the 
velocity and pressure conditions applied on the various boundaries 
occurring In the flow. 

4.1 Free-streaa boundary conditions and Initial values. 

The inflow section In Figure 1, Is Identified with a free-stream 
boundary on the assumption that it Is far enough removed from the bodies 
In the flow so that the Incident uniform flow conditions here are 
unperturbed by the presence of the bodies. Since the free-stream in 
this study is considered to be a uniform parallel flow, constant values 
for velocities and the pressure are always specified on this boundary. 

However, the presence of the accelerating body force g^ in equations 
(3.15) - (3.17) imply that the incident velocities change with time; 
starting from a zero value until they reach their ultimate prescribed 
non-dimensional values. The boundary conditions on F^ are therefore, 


u ■ |g, dt, so that (0 < u 

00 — (XI 

0 

V„ - 0 ; g2 • 0 

«» • 0 J 83 “ ° 

P ■ 0 

00 


< 1) e (0 < t < 1) 


(4.1) 


18 


On the outflow boundaries (r^^, the gradients of the flow vari- 

ables are specified. Denoting the normal to the outflow boundary by 
n, the outflow conditions can be written as 

n • Vu ■ 0 

7P - 0 (4-2) 

4.2 Boundary conditions on body surfaces. 

The flow region Is bounded by a number of body surfaces of 
different families » consisting of the containing walls of the wind 
tunnel and the wing spanning It. The airfoil and the top and the 
bottom walls of the wind tunnel, sections and coincide with 

n ■ constant surfaces- The end wall, k ■ 1 In figure 3, represents 
the C ■ 1 surface. The wall corresponding to C ■ need not be 

considered because the flow Is assumed to be symmetrical about the 
centre-span section at k ■ 4. In the absence of transpiration at these 
material surfaces the velocity boundary conditions belong to the genre 
known as no-sllp boundary conditions for viscous flow. A complete 
absence of relative motion, *slip\ between the solid surfaces and the 
adjoining fluid layers is assumed and the velocity boundary conditions 
are defined by 

u » 0 

v-0 (4.3) 

w * 0 

Boundary values for pressure on body surfaces are calculated 
by designing a momentum equation that is satisfied by the component 
of the pressure gradient normal to the body surface. An expression 
for the normal pressure gradient, which leads to a Neumann type 
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boundary condition, la obtained making uaa of the ralatlona for the 
directional derivative n normal to a surface of constant n and constant 
c (Appendix A) as, 

n • VP - n • (g + ~ V^u) (4.4) 

- ^e ' 

A further simplification results from neglecting the viscous terms in 
the boundary layer and equation (4.4) reduces to 

n • VP “ n • g (4.5) 


The top and the bottom walls of the wind tunnel (P 2 snd F^), are 
coincident with parts of the n ■ boundary, while the airfoil (F^) 

forms a part of the n * surface. These boundaries are therefore 
referred to as n-surfaces. The component of the momentum equation 
normal to a constant n~surface, obtained by replacing n by n^*^^ in 
equation (4.5) is. 


-P Y, + P - - Y 
X 5 y C 


5®1 


(4.6) 


The expression for n^'^\ the normal to constant q-surfaces, is given in 

Appendix A. Using the coordinate transformations for P and P and 

X y 

rearranging the terms in equation (4.6) the Ti“d®rivative of pressure 
on the boundary is obtained as, 

c 

The surface pressure is evaluated from a one sided finite difference 
approximation for 

An expression similar to equation (4,7) is obtained for C-surfaces 
employing the normal to a ^-surface, in equation (4.5). The 

resulting Neumann pressure boundary condition for the k • 1 wall 
(figure 3), is given by 
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(4.8) 



where, 

F'(c) is defined as, 

F'(5) ■ Zj, (4.9) 

The Dlrlchlet velocity boundary conditions (4.1*) together with the 
Neumann pressure boundary condltons (4.7) or (4.8) describe the bound- 
ary values on all material surfaces. 

4.3 Conditions on the re-entrant sections. 

The re-entrant sections, and In figure 1, are not boundaries 
In the physical plane but represent points within the flow field. 

These sections originate from the branch cut from the trailing edge 
of the airfoil to the outflow boundary, made in a C-type grid In order 
to eliminate discontinuity In the geometry of the inner boundary. 
Application of boundary conditions, in the strictest sense, on these 
sections is therefore redundant and not allowed. Maintaining contin- 
uity In the flow variables and their gradients across the cut Is 
de rlgueur. Values of velocity and pressure on these sections should 
therefore emerge as a part of the field solution. However, since the 
re-entrant sections form parts of the n * 1 boundary In the computation- 
al plane, they may be considered subject to a periodic set of boundary 
conditions. 

4.4 Conditions at the trailing edge. 

The trailing edge Is often a sharp point In the physical plane. 
Though the inner boundary, n * 1, is a continuous o.ine In the compu- 
tational plane, the surface-normal vector to the n-surface, n^^^ Is 
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still discontinuous at the trailing edge. The surface pressures at 
the top and the bottom trailing edge points, computed according to 
equation (4.7) are therefore apt to disagree with each other. The 
unequal trailing edge pressures would violate the assumption of absence 
of any unbalanced forces as implied in the no-slip conditions. 

It was deoned advantageous to determine the pressure at the trail- 
ing edge point by means of an average of extrapolates. The extrapolates 
and were calculated from expressions arrived at by representing 
P^ with a three-point one-sided finite difference approximation. Thb 
extrapolates are given by. 


P ■ —(WP - P ) 

ITL 3'' ITL + 1 ITL + 2 ^ 

and 

p m -^4p . P ) 

ITU P ITU-1 ITU - 2^ 


(4.10) 

(4.11) 


ITL and ITU denote the indices in the ^-direction corresponding to the 
lower and the upper trailing edge points respectively. 


4.5 Plane of symmetry conditions. 

As mentioned earlier, computations are performed only on one side 
of the plane of symmetry situated at the mid-span section of the wing. 

In the computational context the second boundary in the z-direction 
needs to be a grid placed beyond the mid-span section so that values 
of the flow variables can be carried over from the grid stationed 
immediately before the mid -span in order to ensure symmetry of the 
flow. Denoting the k-lndex of the centre-span grid by kc, the following 
substitutions are made prior to each SOR iteration. 

(4.12) 


kc + 1 


kc 


\ 
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( 4 . 13 ) 


*’'kc + 1 

V + 1 

^kc + 1 


kc - 1 


■ - w 


kc - 1 



1 


( 4 . 14 ) 

( 4 . 15 ) 
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CHAPTER V 


Fonmilatlon and laplamantatlon of tha Mumarlcal Algor ittaa 

5.1 Obiactlvea and raoulramanta . 

A conciae ravlaw of tha aaavimptiona and objactivaa is parhapt in 
ordar before enbarking upon tha development of the algorithm. The 
governing equations are the Navier-Stokes equations in their non- 
conservative forms. The Poisson pressure equation derives from con- 
siderations of strict conservation of mass; the unconditional satisfac- 
tion of continuity, in effect, furnishing a causal connection linking 
the pressure to the velocities at each temporal state. 

A fully implicit algorithm for solving for the flow variables is 
obtained by means of Backward-Time-Central-Space (BTCS) finite differ- 
encing of derivatives in the transformed plane, thus replacing the 
stability criteria with the less stringent convergence criteria. All 
difference equations apply at the point denoted by the space subscripts 
(l,j,k) and the time superscript (n) . These subscripts and superscripts 
are occasionally omitted for reasons of convenience on the stipulation 
that they are understood to be present in all difference equations. 

The final set of simultaneous algebraic equations is solved by 
a point SOR matrix iterative technique. Special numerical considera- 
lons and computational adjustments will be discussed in their proper 
conte:;ts. 

5 . 2 Coordinate transformations and finite difference approximations . 

The spatial derivatives in the governing equations (3.15 - 3.17), 

in the physical plane are now to be transformed in order to yield the 
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traatfomad equations valid on a rectangular field with sqtiare grid 
in the cMputational plane. For the cake of conciseness of the equa- 
tions, terms such as u , v , etc. will continue to appear in the text 

X y 

but are to be Implicitly assumed to have been evaluated according to 
the relations and definitions given in Appendix A, an well as approxi- 
mated by appropriate finite-difference expressions. The V^( ) terms 
in the momentum and pressure equations are, however, expanded into 
their transformed expressions for a specific purpose to be explained 
in due course. The transformed u-momentum has the following form: 

u^+uu + vu + WU ■-? + - (au_, + 0u + yu.. 

t X y z X ^ j2 ' ^nn CC 

e** 

+ 6u_ + OU, + TU + ^u,) 

5 n ^ t 

+ -^ [2y u + y (u + v ) + u (w + u )1 (5.1) 

R X X y y X *^2 X z'* 


Similar expressions result for the v, w - momentum equations and the 
Poisson pressure equation. 

Space and time derivatives in the transformed Navier -Stokes and 
pressure equations must be approximated by appropriate finite difforence 
quotients on the discrete mesh system. For example, the following 
three-point second-order-accurate central difference approximations Are 
used for subsequent numerical computations: 


3f 

H 


i.J.k 



+ O(AC^) 


(5.2) 



(5.3) 
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with anAlogout fomilAtt for the n and C derivetlvce. It is to be noted 
thet Hi eppeerlng in equetione (5.2) end (5,3) ie chosen to be unity 
by construction. Forwerd end beckwerd three-point finite difference 
expressions were employed to epproxlsiete n*derivetives on the n * 1 nnd 

" \ex bounderies respectively while three-point forwerd expressions 
were used on the C ■ 1 well. These expressions, es well es those for 
Che cross-derlvetives, ere second-order-eccurete expressions in wide 
end conanon use end will not be further eleboreted upon. 

The tine derivetives ere represented by f irst-order-eccurete two- 
point bsckwerd finite difference epproxlnstions et the first step in 
time efter sterting from rest end e second -order , three-point beckwerd 
spproxifflstlon thereefter.' The difference expressions ere given below. 
First-order beckwerds: 

f^” • (-f”"^ + f”)/At (5.4) 

Second -order beckwerds: 

f" - (f**"^ - 4f""^ 4. 3f")/2At (5.5) 

5.3 The Nevier-Stokes equetions in difference form . 

A procedure unique to the prlnltive-verieble fomuletion was used 
to augment the implicit algorithm resulting from the BTCS method of 
finite differencing. The transformed V^( ) terms in the momentum equa- 
tions were expended into their finite difference forms and Che diagonal 
terms, (those with i,Jtk as the spatial subscripts) were separated from 
their off-diagonal counterparts. The diagonal terms were then combined 
with a time-derivative term to form the diagonal elements of the matrix 
representation of the simultaneous algebraic equations. The procedure 
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is best illustrstsd by pttforming the menlpuletion on equetion (5.1). 

The first term in psrentheses on the right hand side of equetion (5.1) 
2 

represents 7 u multiplied by s coefficient. Expansion of the derlvs- 

tivea u.., u and u,, into their finite difference forms result in 

nn c; 

2 

the following expression for 7 u: 

2 1 

’ “ ■ 7 <““« * «“nn * * *“5n + + ♦u^) 

“'“l.J+l.k ■ ^"l,J,k * “l.j-l.k’ 

+ 6u. + CU_ + TU + ♦u,] 

Cn C n ; 


■ 7 -C2. . 26 ^ 2y)Ui ^ kl 

Diagonal 


t 2 ^°^^i-H.J.k ^l-l.J.k^ ^^"i.Jj-l.k ^.J-l,k) 


Off - diagonal 




+ 6u, + OU, + TU + ] 

5n t II ; 


(5.6) 


The tlme-Klerivative on the left hand aide yields one term involving 
on the left hand side and after some rearrangement of terms the 
desired form of the u-momentum equation is obtained as 


n. A ^ 
“ * 


R J 
e 


, (2a + 26 + 2y)] ■ -(uu + vu + wu ) - P„ + 

X y Z A 


R J 
e 


2l»<“l+l,J,k 
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+ 6Up + OU_ + TU + *u,] 

Cn 5 n C 


+ Vx * + V * “I'^x + “t> 1 + h 

e 


(5.7) 


where. 


A - 1 


n - 1 


u 


1 At 


for first-order backwards tlne-diff erenclng. 

(5.8) 


and 


-I 


for second-order backwards time-differencing 
1 ,, n - 1 n - 2, 


_ ^ — A n — ^ V 

‘l ■ 2« “ > 


(5.9) 


V and w-Qomentum equations, developed after a similar fashion, are: 


v”[-^ + — ^ (2a + 26 + 2y)] - -(uv + w -f wv ) 
At jZ X y z 

e 


R j2 ^“^''i+i,j,k Vi,j,k^ ®^''i,j+i,k ''i,j-i,k^ 


^ ^ V ®2 <5.10) 


and, 


^ (2a + 2B + 2y)] * -(uw + vw + ww ) - P 
At R j2 X y z 2 
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+ .iM-i ♦ "i,j ,k-i' * *’'|n * + % + ♦«c> 


+ ^ + u,) + M„(w_ + V.)] + B 


Z Z X 


y y 


(5.11) 


B 2 and B^ are obtained by substituting u with v and w respectively in 
definitions (5.8) and (5.9). 


5.4 The pressure equation. 

The difference form of the Poisson pressure equations is obtained 
by performing the operation described in article 5.3 on the V^P term 
in equation (3.19). The term represents the time-derivative of the 
divergence of the velocity vector. Using a two-point backward differ- 
ence approximation for and noting that d” ■ 0 as dictated by the 
requir^ent for preservation of continuity, the Poisson pressure equa- 
tion is rewritten as. 


2 ^ 2 , 2 , ^ 

« V i I fc,/A u + V + w + 2v u + 2 w u + 2v v 

i»J»K 3C y 2 X y X 2 y 2 


-^(2a + 20 + 2y)F. 




* l'‘^J,k*l ^ ^,3,k-^) * * “‘'l "n * 


2 2 2 2 

- u + u 7 V + p 7 w + p u + p V + p w 

Rg X y ’^z ^xK X yy y ^zz Z 


+ p (u +v)+p (w +u)+p (v + w)l- 
^xy y x' *^zx x z' ^zy' z y'-' At 


D° ' ^ (5.12) 


5.5 Difference equations for boundary conditions. 

Pressure values on impermeable surfaces In the flow are determined 
by extrapolation from the field subject to the Neumann boundary condi- 
tions delineated in article 4.2. P or P in equations (4.7 - 4.8) are 

n ^ 
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evaluated on the surfaces by three»polnt backward or forward difference 
approximations as the case may be, depending on the location of the 
surface with respect to the coordinate system. 

The n ~ 1 boundary is congruent in part with the airfoil; there- 
fore ->-derlvatives on the airfoil are of necessity to be approximated 
by forward differences. Along a similar line of reasoning, the top 
and the bottom walls (n ■ \ax^ ' seen to require backward n- 
differencing while the end-%rall = 1) requires forward ^-differencing. 

The resulting expressions for boundary pressure on the various 
solid boundaries are: 

On the airfoil, 

1 4 2 1 

i,l,k r i,3,k ^ r i,2,k 3 Yj. 5 c 

+ (5.13) 


On the top and the bottom walls. 


^i,JMAX,k 


- ±P + ^ 2 1 

3 i,JMAX-2,k 3 i,JMAX-l,k + ^ — [P^S 

-> Yc ^ c 

+ J,.(g2*5 ■ (5.14) 


On the end-wall ( ? = 1) , 

"^i,j,i yi,j,3 ^ri,j,2 


F (5)g3 


(5.15) 


5.6 Special consideritions and computational adjustments. 

5.6.1 Finite differencing on the reentrant sections. 

The procedure for the application of finite difference approxima- 
tions on the cut extending from the trailing edge of the airfoil to 
the outflow section in a C-type grid is rather extraordinary and 
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deserves special attention. The t\to reentrant sections, and F^ 

(fig. 1), resulting from the cut are one and the same line in the phy« 
sical plane, albeit ostensibly separated from each other in the trans- 
formed computational plane. Any point on the segment F^ is congruent 
with a corresponding point on F^, eg., the points designated by the 
indices (IMAX-1, 1) and(ITV + 1, 1) are indistinguishable in the phy- 
sical plane from the points at (2,1) and (XXL - 1, 1} respectively with 
similar correspondence existing between all other points distributed 
between the trailing edge and the outflow boundary. These correspond- 
ing points, therefore, have the same x and y coordinate values inspite 
of the differences in their associated C-values. Evaluation of deriv- 
atives and flow calculations are necessary along only one of these 
reentrant sections since they carry the same values of the flow vari- 
ables. 

At a casual glarxe, the evaluation of derivatives on either of 
these sections by means of central-difference approximations may se^ 
rendered impossible by the absence of a J - 1 line, since the reentrant 
sections correspond to the n * 1 line. A closer view of the physical 
geometry involved reveals, however, that the J * 2 line reemerging 
behind the trailing edge after circumscribing the airfoil, serves as 
a surrogate J - 1 line thus making central-differencing possible while 
ensuring satisfaction of the strict requirement for continuous deriva- 
tives across the cut. The resulting central-difference expressions 
for the derivatives across the cut at the points corresponding to 
I • II and 1 = 12 (figure 1) are given below. 



( 5 . 17 ) 


^^nn^il.l “ ^^nn^i2,l “ ^11,2 “ ^^il,l ■•■ ^ 12 , 2 

^^e.n^ii.i “ ^^ 5 n^i 2 ,i = ^^ 11 - 1,2 ■ ^ 11 + 1,2 ■*■ ^ 12 + 1,2 

■ ^ 12 - 1 , 2)/A (5.18) 

where » 

12 * IMAX - (II - 1) 

The k-index has been omitted for reasons for simplicity* 

5.6.2 Artifical viscosity. 

Extraneous oscillations contaminating the solution, eg., those 
arising out of inaccuracies introduced by finite differencing, may 
lead to instability and ultimately divergence of the solution particu- 
larly at high Reynolds numbers. An artificial viscosity, judiciously 
used, can diffuse these nonphysical instabilities and save the solution 
from being overcome by violent oscillations. 

Since the preservation of continuity is an underlying stipulation 
in the solution scheme, an unusually high magnitude of the divergence 
V • u of the velocity vector should be symptomatic of spurious oscill- 
ations in the dependent variables. Based on this reasoning an artifi- 
cial viscosity coefficient, [14], designed to be proportional to the magni- 
tude of V • u, is appended to the viscosity coefficient in the momentum 
equations (3.15 - 3.17). The coefficient p is therefore composed of 
three parts as shown below, 

U » 1 + e + d (5.19) 

where, 

c accounts for the turbulent eddy viscosity 
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and 

6 represents the artificial viscosity being discussed here 
6 is being defined as 

0 - XR jIv • ul (5.20) 

e ' 

where X ^ 0 and J is the Jacobian of the coordinate transformation. 

The artificial viscosity 6 can be perceived to have the dual 
sough'.-after attributes of coming into play only where 7 • u is un- 
desirably high and of being proportional to the magnitude |v • u| at 
the same time. Suitably chosen values of 1 can limit the effects of 
the artificial viscosity. 


5.7 A brief outline of computational sequences. 

The wing, the wind-tunnel and the associated coordinate system are 
started from rest and accelerated to the ultimate free-stream velocity. 
The free-stream velocities are distributed in time along a cosine curve 
of the following form. 


U ■ 4[1 + cosir (t + 1)] 0<t<l 

00 / — 

U - 1 for t > 1 

CO — 


(5.21) 


At each time-step the boundary values of pressure and the veloci- 
ties are applied to the various boundaries occurring in the flow prior 
to each iteration. One iteration of the SOR iterative method is then 
performed. Appropriate optimum acceleration parameters co, are calcu- 
lated as described in Appendix B, and used to update the successive 
Iterates according to the following formula. 


.(s + 1) 


f'“ • - 0)f (1 -0))f^®^ (5.22) 

where f represents any one of the three dependent variables u, v and 
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(It, lAile th« wpcracrlpt (•) denotas the itcratlona counter. The 
preeeure eolutlon le not accelerated. The unaccelerated Gauas-Seldel 
Iterate is represented by 7. 

The Navler-Stokes and the pressure equations are solved sijnultan- 
eously In the field to maintain the implicit characteristic of the 
algorithm. Iteration and evaluation of boundary conditions are alter- 
nated until a desired degree of convergence is achieved at the current 
time level. The free-stream is then accelerated to the next step in 
time, the boundary /alues are calculated and a new set of iterations 
commence. This procedure is continued until a steady state solution 
results. 
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CHAPTEK VI 


CoBputatloiutl Results 
6.1 Coordinate Systems 

Coordinate grids were generated £or symmetric NACA airfoils for 
a range of Reynolds numbers by means of the boundary-fitted schemes 
of Thompson [2] and Sorenson [3]. The airfoils used were the NACA 
0012 and NACA 66^018 sections. Field slses of the computational grids 
varied with variations In the Reynolds number and the shape of the 
airfoil used. 

In the Reynolds number range of 1000-10,000, each two-dimensional 
grid contained 49 x 22 points describing the field geometry, with five 
such grids placed along the span at stations Including the side-wall 
and the post-mldspan section. In other words, there were 49 constant- 
C lines and 22 constant-n lines defining each discrete two dimensional 
mesh. 

For R » 40.000. the two dimensional field size was Increased to 
e 

55 X 31 to accommodate the intensified requirements for resolution of 
the boundary layers. The number of spanvise stations was increased to 
20 keeping in view a proposed project involving the resolution of the 
boundary layer near the side-wall at high Reynolds numbers. 

Coordinate lines were concentrated in the boundary layers, the 
first line away from the surface being located at a distance of one 
percent of the calculated thickness of the boundary layer on a flat 
plate. 

The computational results presented in this dissertaion were ob- 
tained from computations performed on grids generated by Sorenson^ s 
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'GRAPE' code, Although grid generetloc teehniqiiee of both T b oa p e ott and 
SereMcm were eatpegfae nt e d with. Picwroe 4*7 fXMent the coordinate 
eyitems for the particular flow situations Investigated. In all these 
coordinate grids the chord of the airfoil lies midway between the top 
and the bottom walls of the wind tunnel along the y ■ 0 line with the 
leading edge at x ■ 0, and the trailing edge at x ■ 1. The top and 
the bott<»n walls correspond to y * 2, and y >2 lines and have their 
leading emd trailing edges at x ■ -1 *0, and x ■ 4 • 0 respectively. 

The span extends from s - 0 to z - 2, so that the side-wall corresponda 
to the z ■ 0 plane while the mldspan section Is situated at z > 1. 

6.2 Results for NACA 0012 Airfoil 

6.2.1 Results for R ■ 1000 and A ■ 0* 

e 

Relevant flow condltons governing this part of the computations 

are: 

K - 1000 
e 

A - 0* 

0 ■ 0; X • 0 
e : Computed 

The turbulent eddy viscosity e, though inappreciably small in 
magnitude at this low Reynolds number, was included in these computa- 
tions but was quiescent as might be expected. The use of the artificial 
viscosity 6 was neither warranted nor necessary since no instability 
was encountered in the solution. 6 was set equal to zero by using a 
zero value for the coefficient X. 

Distributions of C^, the pressure coefficient are presented at 
two time levels - t - 1.0 and t <■ 2.0. The transient states of the 
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solutions, though unimportant from tha point of vlaw of practical 
utility, ara prasantad to daaonstrata tha tima-dapandanca of tha sol- 
ution. Figures 8-11 show the C distributions at t ■ 1.0, i.a., after 

P 

the solution has advanced through 100 time steps, each .01 in magni- 
tude. C distributions In figures 12-15 are those at t ■ 2.0. A 
P 

distinct change in the pressure profiles is noticed particularly near 
the leading and the trailing edges. The differences between the 
pressure values on tha upper and the lower surfaces are seen to have 
considerably lessened. These discrepancies arise out of inadequate 
resolution of the flow quantities at regions containing large gradients 
and are attributable in part to an insufficient number of points in 
the discrete mesh. 

Velocity prof lies around the airfoil are represented in figures 
16-19, while figures 20-23 show the distribution of velocities along 
constant-^ lines bridging the gap between the airfoil and the top wall 
of the wind tunnel. Well developed boundary layers are noticed at both 
surfaces. 

Resolution of the boundary layer adjacent to the side-wall (C “ 1) 
was not attempted, the number of spanwlse grids being limited by con- 
siderations of affordability of computer storage. However, pressure 
profiles do vary along the span as is evident upon examination of the 

C distributions at different values of the index k. Little or no 
P 

change was noticed in the pressure distributions on the top and the 
bottom walls with change in the value of k, and the Cp-^iistributlons 
at k ■ 4 are presented In figures 24 and 25. 

The solution reached a reasonable degree of convergence and was 
considered to have demonstrated a proper qualitative trend at t ■ 2.0 
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where the celculetlone were termlneted. 

6.2.2 Reeulte for ■ 1000 end A ■ 12* 

The coordlnete syetem shown In Figure 5 wee used In this pert of 
the computet ions. The relevent flow peremeters were, 

R^ - 1000 

A - 12* 

0 - 0; X - 0 

e : computed 

A greet deel of computet ion-time wee seved In this cese by sterting 
the solution from the steady-stete solution obtained in the case of 
zero-angle-of-attack described in the previous section. 

Velocity profiles at the midspan section (k ■ A) are presented 
in Figures 26 and 27 at 5 and 60 time-steps after restart respectively. 
The beginnings of a vortex, arising out of the sudden rotation of the 
airfoil, is discernible at the trailing edge; nevertheless, the solu- 
tion scheme demonstrates remarkable adaptability and the flow soon 
settles down. The velocities follow the contours properly and a stable 
solution is obtained as shown in Figure 27. 

Pressure coefficient distributions along the span are presented 
in Figures 28-31 and generally agreeable qualitative trends are observed. 

6.3 Results for NACA 66^018 airfoil section. 

6.3.1 Results for R, • 10,000 and A • 0* 

a [ . 

The governing parameters are: 

Rg - 10,000 

A ■ 0 
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X - 1. 0 


e : computad 

This phase of the computation was performed on the mesh shown in 
Figure 6. The solution was restarted using the steady state eolution 
obtained with the NACA 0012 airfoil described In section 6.2.1, as 
an inital guess. 

As the solution progressed in time the pressure went through 
phases of develoiH&ent and finally settled to a steady profile after 
90 time steps. Figures 32-3A pres«it the temporal evolution of the 
pressure distribution on the airfoil at the k ■ 4, station, as the 
flow strive to onform to the suddenly imposed changes In geometry 
and Reynolds number. Comparing the results with those obtained in 
the case of the NACA 0012 airfoil at A * 12”, the change in the 
Reynolds number seemed to be the harder obstacle to overcome. The 
artificial viscosity was brought into effect by setting 1. ~ 

distributions at various spanwise stations on the airfoil are presented 
in figures 35-38. 

Experimental results obtained by Mueller [l5l, from wind-tunnel 
tests of the NACA 66^018 airfoil at a Reynolds number of 40,000 are 
plotted along with the computational results for R^ * 10,000, in Figure 
39. Since the Reynolds numbers of the results being compared are 
different, no degree of quantitative exactitude of replicability was 
sought. The comparison was performed in the hope of establishing that 
the computed pressures do duplicate the general qualitative trends df 
the experimental results and this assertion seems to have been borne 
out by the plotted comparison. 
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Th« lack of a subatantlal dagraa of staepnaas in tha praaaura 
proflla at tha laading adga la attributabla in part to Inauff Iclancy 
of tha nunbar of grid polnta avallabla in thla raglon of aharp grad- 
lanta in flow variablaa. Tha trailing adga praaaura, howaivar, dropa 
with tlna and Indlcataa a promlaa of following tha trand aat by tha 
exparlnantal raaulta. 

Tha velocity profilaa on tha aurface of tha airfoil at tha aid- 
apan aaction ara praaantad in Figure 40. Tha boundary layer la aaen 
to be wall devalopad. Tha apparent penetration of tha noaa of tha 
airfoil by tha fluid reaultad from the abrupt changaa in geometry and 
Rejmolda number compounded with altered grid apacing in tha boundary 
layer as well aa the inadequacy of tha number of n-linaa available 
near the aurface in properly resolving the boundary layer. Continued 
advance in time did not get rid of this abnormality. 



CHAPTER VII 


Dlscuaslons and Concluaions. 

Tha lapllcit algoritha anployad in this raaearch to computa tha 
flow about an airfoil spanning a wind tunnal doas liva up to most of 
tha raaaarchar'a axpactations although tha absanca of quantitative 
comparison of results leave some things to be desired. 

A scrutiny of tha presented computational results Indicates that 
an algoritha with substantially stable characteristics is at work. The 
stability, expected of a fully implicit algorlthir., asserts itself, 
particularly in the case of the NACA 0012 airfoil at an angle of attack 
of 12 degrees. Despite restarting the solution from the steady state 
obtained with zero angle of attack, the flow soon attains conformity 
with the changed geometric configuration. The front stagnation point 
moves away from the leading edge to the adjacent point cn the lower 
surface of the airfoil, the streamlines follow the contours and the 
pressure builds up on the upper surface while dropping on Che lower 
surface. The flow demonscrates similar, if not as spectacular, adapt* 
ability in the case of the NACA 66^018 airfoil with “ 10,000. The 
results, therefore, are encouraging enough Co warrant planning Che use 
of the algorithm in future cominitaClons at higher Reynolds numbera. 

Proper resolution of a flow field is directly dependent on the 
number of grid points available in the field. It was in this crucial 
area that the research was constrained by the limits of affordability 
of computer storage and hence the number of discrete mesh points. The 
storage requltements and the computation time for a complete three 
dimensional solution are truly overwhelming and consequently it was 
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deemed prudent to limit the field size to one Just large enough to make 
It possible to ascertain the qualitative characteristics of the solu- 
tion. 

Some of Che relevant computational parameters should be made note 
of here. The number of iterations for each time-step was limited to a 
mflylimim of 20. Further iterations did not seem to significantly en- 
hance convergence. Typical error norms after 20 iterations were of the 
order of 10 • With a total of 4y x 22 x 5 points in the three-dijnen- 

sional computational space, an average of 50,0 seconds of computer time 
per time^step was expended on the GDC CYBER 175 machine. The storage 
requirements amounted to an approximate total of 3 x lOg locations. 

With a limit of lOOK for the machine used, there was not sufficient 
free space left to enable an increase of the number of data points. 

Among the forecasts for the future is a project already in progress, 
involving computations at realistically high Re3molds numbers, on a 
grid many times bigger and denser than the ones used in this research. 
The GDC GYBER 203 offers virtually unlimited storage capability and is 
the machine chosen for these future computations. 
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Figure 4. Coordinate System for NACA 0012 Airfoil - 

R = 1000, A = 0“ 
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Figure 5. Coordinate System for NACA 0012 Airfoil 

R = 1000, A = 12“ 
e 
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Figure 6. Coordinate System for NACA 66j018 Airfoil — 

R = 10,000, A * 0“ 
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Figure 7. Coordinate System for NACA 66j018 Airfoil - 
= 40,000, A - 0“ 
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Figure 8, Pressure Distribution (NACA 0012) - 

R « 1000, A * 0% t - 1.0, k « 1 (wall) 
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Figure 9. Pressure Distribution (NACA 0012) - 

R = 1000, A = 0°, t = 1.0, k » 2 
e 
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Figure 10. Pressure Distribtuion (NACA 0012) 
Rg - 1000, A * 0", t - 1.0, k - 3 



-I4>&iiii!niliiijijiiiitiiiniiilit inini!iimiiii i 

0 •! -6 t 

X 


Figure 11. Pressure Distribution (NACA 0012) 

R = 1000, A = 0% t » 1.0, k = 4 

e ’ 


(tnidspan) 


51 



ORIGINAL PAGE IS 
OF POOR QUALITY 



Figure 12. Pressure Distribution (NACA 0012) 
Rg - 1000, A - 0*, t - 2.0, k - 1 



Figure 13. Pressure Distribution (NACA 0012) 
- 1000. A - 0*, t - 2.0, k - 2 


(wall) 
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Figure 14. Pressure Distribution (NACA 0012) - 
- 1000, A » 0% t - 2.0, k - 3 
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Figure 15. Pressure Distribution (NACA 0012) - 

» 1000, A » 0®, t » 2.0, k * 4 (midspan) 
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Figure 16. Velocity Profiles Around Airfoil (NACA 0012) 
- 1000, A - 0", t • 2.0, k - 2 



Figure 17. Velocity Profiles Around Airfoil (NACA 0012) - 

R • 1000, A • 0“, t - 2.0, k - 3 
e 
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Figure 18. Velocity Profiles Around Airfoil (NACA 001?) 

■ 1000, A«0®, t»2.0, k*4 (midspan) 



Figure 19. Velocity Profiles Around Airfoil (NACA 0012) 
- 1000, A - 0®, t - 2.0, k - 5 
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Figure 22. Velocity Distribution in Wind Tunnel (NACA 0012) 
Rg - 1000, A = 0% t = 2.0, k » 4 (midspan) 



Figure 23. Velocity Distribution in Wind Tunnel (NACA 0021) 
Rg = 1000, A = 0% t = 2.0, k = 5 
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Figure 24. Pressure Distribution On The Top Wall - 

R = 1000 
e 
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Figure 25. Pressure Distribution On The Bottom Wall - 

R = 1000 
e 
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Figure 26. Velocity Profiles Around Airfoil (NACA 0012) 

R = 1000, A = 12% t = 2.05, k = 4 (midspan) 
e 
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Rg = 1000, A = 12“, t - 2.6, k = A (midspan) 
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Fieure 28. Pressure Distribution (NACA 0012) - 

R = 1000, A = 12°, t = 2.05, k = 4 (midspan) 
e 



Figure 29. Pressure Distribution (NACA 0012) - 

R = 1000, A= 12° t = 2.1, k= 4 (midspan) 


0 


.4 


X 


3 


1-0 


Figure 31. Pressure Distribution (NACA 0012) - 

R = 1000, A = 12®, t = 2.6, k = 4 (midspan) 
e 
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Figure 32. 




Figure 34. 


Pressure Distribution 
(NACA 663018 ) 

- 10»000, A • 0* 

t - 2.1, k - 4 


Pressure Distribution 
(NACA 663OI8) 

Rg - 10,000, A - 0* 

t - 2.2, k - 4 


Pressure Distribution 
(NACA 66^018) - 
Rg = 10,000, A = 0“ 

t = 2.3, k = 4 
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Figure 35. Pressure Distribution (NACA 662 OI 8 ) 

= 10,000, A « 0", t = 2.9, k * 1 (wall) 
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Fiiure 37. Pressure Distribution (NACA 663 CI 8 ) - 
R - 10,000, A * 0®, t - 2.9, K • J 
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Figure 38. Pressure Distribution (NACA 665 OI 8 ) - 

Rg * 10,000, A*0", t*2.9, k»4 (tnidspan) 
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Figure 39. 


Comparison of Computational Results 6^3018) 

at R - 10,000 with Experimental Results a 


'40.000, A 


0 ® 
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Figure 40. Velocity Profiles (NACA 66^018) - 

- 10,000, A - 0*, t - 2.9, k - 4 (midspan) 
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APPENDIX A 


Relations and definitions in the transformed plane. 

Contained herein are the pertinent relationships defining the 
transformations of spatial derivatives, appearing in the governing 
equations, from the physical to the computational space. All the 
transformations are given in their fully non-conservative forms. The 
following notations for scalar and vector functions are valid through- 
out all the expressions. 

f(x,y,z,t) - A scalar funclton with at least two continuous 
derivatives. 


F(x,y,z) = iFj^(x,y,z) + jF2(x,y,z) + k(x,y,z) - A continuously 


differentiable vector-valued function. j» and 

k denote the Cartesian unit vectors. 

Definitions of the transformations in 2-D grid generation. 

2 ^ 2 
a = X + y 

(A.l) 

c n n 

= x.x + y^y 

c n n 

(A.2) 

2 ^ 2 

(A.3) 

■"c = 

(A.4) 

Transformations in the 

three-dimensional space. 

Z = F(;) 

(A. 5) 

'2 

a = a^F {0 

(A. 6) 
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»2 


8 • y^f (0 

(A.7) 

CM 

O 

*n 

N 

>- 

(A.8) 

.2 

(A. 9) 

t 

” ■ j‘' ‘%'v - %“*> 

(A. 10) 

t 

T “ (y^i>x - x^Dy) 

(A. 11) 

II 

♦ ■ - Vs' 

(A.12) 

J - 2^(xjy„ - yjx^) - f’u).^ 

(A. 13) 

• “sc * ®*6n * *”nS 

(A. 14) 

Dy - ayj; + + 6y^^ 

(A. 15) 


Der ivat ive transformations 




(A. 16) 


ty - - <Vn • 


(A.17) 


f - (3f/3z)^ ^ ^ - f(;)/F iO 


(A. 18) 


' W«a ^ h\M 

2 '> 3 

* hr * h^m'> <%'« - 

* ^>’c='r,’‘Cn * ‘ <*•!’> 
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f 

yy 



2 


* ■ ^5%”?-. * ‘«Sn>‘’n'c 




2 2 
+ (x x._ - 2Xj.x X + X X„ )(y-f_ 

n 55 5 n 5n 5 nn 5 n 

- Ve>'^c’ 

(A. 20) 

f 

= 0^f/3z^) . = [f__ - f_F"(5)/F'(?)l/F'^(C) 

(A. 21) 

zz 

x,y,t 



f 

xy 

• 'Vn * ’‘r?(>Ur, - *?>'c'nn ‘ 



* - ‘*5^ * * W 


^5%>V 


* ■ <’‘5% * *n>’5>»«n * *C>' 

C^r,' 



(x,f - X f )/J ^ 

5 n n 5 c 


(A.22) 

f 

= (f X^ - f^-X )/J 


(A.23) 

yz 

n; 5 55 n 



f 

= (f^ V - f vj/j 


(A.24) 

xz 

CC-'n n5'^5 




The Laplacian Operator. 




■ ‘“'55 «*nn * ‘‘sn “'c ^ 

xfn 





(A.25) 


Unit Normal Vectors. 




Nonoal to ri-surface: 

i Vn/jvnl * (“iy^ + (/“i. 26 ) 

Normal to 5-surface: 

i V 5 /I 75 I = (iy^ - (A. 27) 
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Directional Derlvat Ives . 


3£/3n^'^^ = 

• V ■ ‘I'c'n ■ 

(A. 28) 


• !' ■ <“c^ - 

(A. 29) 

3f/3n^^> E 

• 7f - f./F’(0 

(A. 30) 
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APPENDIX B 


Optlnuo acceleration parameters. 

The Navier-Stokes equations In non-conservative primitive varia- 
ble formulation fit the description of the following general partial 
differential equation In the transformed plane. 

*1^5 * * “I'c * ‘ 2 'n * 's's + <»•» 

f identifies itself with u, v or w in the x, y or znaomentum equation* 
The spectral radius for Jacobi iteraion of a set of IMAX x JMAX 
X KMAX simulataneous difference equations with constant coefficients [14] 
is given by, 

|C - 2(A^ + A 2 + A^) I |4A^ - I ^ 

- B.^l Cos( )] 

KMAX + 1 

(B.2) 

for SOR Iteration are easily 
expression, 

2 

4A2^ > ^^ 3 ^ - ®3 

(B.4) 

J 

The coefficients in equations (B.l) are listed on the following page. 


CosC jtJ ^ 


The optimum acceleration parameters w, 
obtained utilizing in the following 


2 2 2 
if4A^^>B^^, 




and, 


otherwise. 


1 + "TTr^ 
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Coefficients in the x-momentum equation: 



(B.5) 

Aj - Rj 

(B. 6 ) 

*3 • *3 

(B.7) 

“l ■ ^1 + E j Vn ■ Vn* 

e c ^ ' 

(B. 8 ) 

*2 ■ ■^2 + R J <- Vs’ 

e c 

(B.9) 

1 

■ ■'3 ^ K .V . 

e F U) 

(B.IO) 

T 

c * - ^ 

At 

(B.ll) 

Coefficients in the y-momentum equation: 



(B.12) 

A 2 .R 2 

(B.13) 

Aj - R 3 

(B.14) 

“1 ■ ■'V R^ <- ^Vn * "xV 

(B.15) 

‘2 ■ ^2 ♦ - ^ 25 ) 

(B.16) 

^ ^ F (;) 

(B.17) 

T 

G * - 

At 

(B.18) 


Coefficients in the z-momentum equation: 

\ \ (B.19) 
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Aj - Rj 
A3 - R3 


'1 ■ ■'1 lA: ‘ Vn ■ 


e c 


“2 ‘ ■'2 * 57; <- Vs * Vc* 


2u 

®3-''3 + ir 


z 1 


e F (O 




The repeated terms in the above expressions are; 


V 


R, 


2 2 

V 


% 


R J 
e 


T - uy - vx ) + — ^ a 

1 Jc n ^ R j2 




T3 - -r^ + * 

^ F (;) R^r 


T ■ 1 for first-order time differencing 


T ■ for second-order time differencing 


(B.20) 

(B.21) 

(B.22) 

(B.23) 

(B.24) 

(B.25) 

(B.26) 

(B.27) 

(B.28) 

(B.29) 

(B.30) 

(B.31) 

(B.32) 

(B.33) 
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